home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / poly.r < prev    next >
Encoding:
Text File  |  1994-09-20  |  1.7 KB  |  72 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2.  
  3. //  Syntax: poly(A)
  4.  
  5. //  Description:
  6.  
  7. //  Find characteristic polynomial of A
  8. //
  9. //  If A is an n by n matrix, poly(A) is an n+1 element row vector
  10. //  whose elements are the coefficients of the characteristic polynomial
  11. //      det(A-zI) = 0.
  12. //  The output coefficients are stored in decending powers as
  13. //      [C1,C2,....,Cn+1]
  14. //  It represents polynomial
  15. //      C[1]*z^n + C[2]*z^(n-1) + .... + C[n]*z + C[n+1]  
  16. //
  17. //  If A is an n by 1 column vector, elements of A are roots of 
  18. //  polynomial
  19. //      (z-A[1])(z-A[2]) ... (z-A[n]) = 0
  20. //  The output coefficients represents the expanded form of above
  21. //  polynomial.
  22.  
  23. //  Example:
  24. //  > X=[5,-6;1,0]
  25. //   X =
  26. //          5         -6
  27. //          1          0
  28. //  > p=poly(X)
  29. //   p =
  30. //          1         -5          6
  31. //  > r=roots(p)
  32. //   r =
  33. //          2
  34. //          3
  35. //  > p1=poly(r)
  36. //   p1 =
  37. //          1         -5          6
  38. //
  39.  
  40. //  See Also: roots
  41. //
  42. //  Tzong-Shuoh Yang (tsyang@ce.berkeley.edu)  5/7/94
  43. //
  44. //-------------------------------------------------------------------//
  45.  
  46. poly = function(A)
  47. { 
  48.     if (A.nc == A.nr) {
  49.       // A is a matrix, find its eigenvalues
  50.       e = eig(A).val';
  51.       l = inf();
  52.       // trim INFs
  53.       e = e[find(e != l && e != -l)];
  54.     else if (A.nr == 1 || A.nc == 1) {
  55.       // elements of A are roots of the polynomial
  56.       e = A[:];    // force a column vector
  57.     else
  58.       error("Argument of poly() must be a square matrix or a vector");
  59.     }}
  60.     // construct char. polynomial
  61.     c = [1,zeros(1,e.nr)];
  62.     for (i in 1:e.nr) {
  63.        j = i + 1;
  64.        c[2:j] = c[2:j] - e[i]*c[1:i];
  65.     }
  66.         
  67.     if (all(imag(c) == 0)) {
  68.        c = real(c);
  69.     }
  70.     return c;
  71. };
  72.